%Create 
FileNamesOuter=cell(1,6);
FileNamesOuter{1,1} = {'FirstWorkspace6AH20','SecondWorkspace6AM20','ThirdWorkspace6AU20','FourthWorkspace6AZ20','FifthWorkspace6AH21'};
FileNamesOuter{1,2} = {'FirstWorkspace6BH20','SecondWorkspace6BM20','ThirdWorkspace6BU20','FourthWorkspace6BZ20','FifthWorkspace6BH21'};
FileNamesOuter{1,3} = {'FirstWorkspace6CH20','SecondWorkspace6CM20','ThirdWorkspace6CU20','FourthWorkspace6CZ20','FifthWorkspace6CH21'};
FileNamesOuter{1,4} = {'FirstWorkspace6EH20','SecondWorkspace6EM20','ThirdWorkspace6EU20','FourthWorkspace6EZ20','FifthWorkspace6EH21'};
FileNamesOuter{1,5} = {'FirstWorkspace6JH20','SecondWorkspace6JM20','ThirdWorkspace6JU20','FourthWorkspace6JZ20','FifthWorkspace6JH21'};    
FileNamesOuter{1,6} = {'FirstWorkspace6SH20','SecondWorkspace6SM20','ThirdWorkspace6SU20','FourthWorkspace6SZ20','FifthWorkspace6SH21'};
%Create string with root letters for different markets, A=AUD/USD,
%B=GBP/USD, C=CAD/USD, E=EUR/USD, J=JPY/USD,S=CHF/USD
Strcatroots = {'A','B','C','E','J','S'};

%Change directories
functionsFolder=cd('..');
parentDirectory=cd('output\');


nVec = [1];
outerloop = 0;
%outerloop cycles through different roots (contract markets)
%can change outerloop to 6 from 5 to run CHF/USD test
while(outerloop<5)
    tic
    outerloop = outerloop + 1;
    FileNames = FileNamesOuter{1,outerloop};
    n = nVec(1);

    %Initialize values for each storage variable
    
    %Beta Coefficients for different models and standard errors
    BetasMech = cell(11,1);
    BetasMech05d = cell(11,1);
    BetasNMech = cell(11,1);
    StandardErrorsMech = cell(11,1);
    StandardErrorsMech05d = cell(11,1);
    StandardErrorsNMech = cell(11,1);
    
    %Define general nonlinear model to estimate
    nlModel = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^Beta(3)).*sign(X))); 
    %Define square root price impact model
    nlModel2 = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^0.5).*sign(X))); 
    

    
    LambdasNMechL = cell(11,1);
    BetasNMechL = cell(11,1);
    SandwichesNMechL = cell(11,1);
    DaySizesNMechL = cell(11,1);
    R2sNMechL = cell(11,1);

    LambdasMechL = cell(11,1);
    BetasMechL = cell(11,1);
    SandwichesMechL = cell(11,1);
    DaySizesMechL = cell(11,1);
    R2sMechL = cell(11,1);
    VolumesMechL = cell(11,1);
    AdditionalValues = cell(11,1);

    DatesForLater = cell(11,1);

    %Four vectors to build out: price with mechanical depletion, price with
    %only non-mechanical, signed order size, date/time (Y1Mech Y1NMech X1 t2)
    Y1MechAll = [];
    Y1NMechAll = [];
    X1All = [];
    t2All = [];
    retsAll = [];

    filecounter = 0;
    %Cycle through each contract month per market: Mar '20, June '20, Sep
    %'20, Dec '20, Mar '21
    while(filecounter<5)
        filecounter = filecounter + 1;
        clear ESU1820180801;
        
        %Load table of data for contract month
        load(FileNames{1,filecounter});
        
        %Convert timezome from UTC to NYC time
        a2 = datetime(ESU1820180801{:,1},'TimeZone','UTC');
        a2.TimeZone = 'America/New_York';
        t2 = datenum(a2);
        
        %Clean up NaNs if any
        dataIsNaNFormat = double(ESU1820180801{:,2:7});
        dataIsNaNFormat(dataIsNaNFormat == 0)= NaN;
        
        %Run routine that analyzes dataset to create a list of aggressive
        %orders signed by direction, with size and order book dynamics
        %corresponding to order time
        [AggressiveOrders2, ~,Y2] = MicroStructureES([t2 dataIsNaNFormat]);
        
        %Add a column that consists of the "true price"
        AggressiveOrders2(:,13) = (AggressiveOrders2(:,8).*AggressiveOrders2(:,9) + AggressiveOrders2(:,10).*AggressiveOrders2(:,7))./(AggressiveOrders2(:,8) + AggressiveOrders2(:,10));
        
        %Filter out trades outside of liquid trading hours
        t2H = hour(AggressiveOrders2(:,1));
        t2M = minute(AggressiveOrders2(:,1));
        
        selector = ((t2H>=18) +(t2H<15)); %FX & US Rates
        %selector = ((t2H>=3) +(t2H<12)); %Europe Rates
        AggressiveOrders3 = AggressiveOrders2(selector>0.01,:);
        clear AggressiveOrders2;
        
        %Create X and Y vectors using careful elimination of intraday jumps
        t2H = hour(AggressiveOrders3(:,1));
        diffjumps = ((abs(t2H(2:end)-t2H(1:end-1))>2)+ (abs(t2H(2:end)-t2H(1:end-1))<22))>1.01;
        diffjumps = [diffjumps; 0]; %#ok<AGROW> %next order is following day
        
        %Determine number of days in the dataset
        selectorDays = diffjumps.*(1:size(diffjumps,1))';
        selectorDays = [0;selectorDays(selectorDays>0,:);size(AggressiveOrders3,1)];
        
        %Pre-allocate based on number of days for daily estimates
        betasMech = zeros(size(selectorDays,1)-1,3);
        betasMech05d = zeros(size(selectorDays,1)-1,2);
        betasNMech = zeros(size(selectorDays,1)-1,3);
        standardErrorsMech = zeros(size(selectorDays,1)-1,3);
        standardErrorsMech05d = zeros(size(selectorDays,1)-1,2);
        standardErrorsNMech = zeros(size(selectorDays,1)-1,3);
        lambdasNMech = zeros(size(selectorDays,1)-1,1);
        betasNMechL = zeros(size(selectorDays,1)-1,2);
        sandwichesNMech = zeros(size(selectorDays,1)-1,4);
        daysizesNMech = zeros(size(selectorDays,1)-1,1);
        r2sNMech = zeros(size(selectorDays,1)-1,1);

        lambdasMechL = zeros(size(selectorDays,1)-1,1);
        betasMechL = zeros(size(selectorDays,1)-1,2);
        sandwichesMechL = zeros(size(selectorDays,1)-1,4);
        daysizesMechL = zeros(size(selectorDays,1)-1,1);
        r2sMechL = zeros(size(selectorDays,1)-1,1);
        volumesMechL = zeros(size(selectorDays,1)-1,1);
        additionalValues = zeros(size(selectorDays,1)-1,6);
        
        datesForLater = zeros(size(selectorDays,1)-1,1);
        
        counter = 0;
        while(counter<size(selectorDays,1)-1)
            %Loop through days and run analysis
            counter = counter + 1;
            startval = selectorDays(counter,1)+1;
            endval = selectorDays(counter+1,1);
            datesForLater(counter,1) = AggressiveOrders3(startval+1,1);
            %n=1; %n controls the number of transactions ahead for which we
            %are calculating price impact
            %X1 is signed order size
            X1 = AggressiveOrders3(startval:endval-n,3).*AggressiveOrders3(startval:endval-n,4);
            %Y1Mech is the difference in true price from before one trade
            %to the subsequent trade
            Y1Mech = AggressiveOrders3(startval+n:endval,6)-AggressiveOrders3(startval:endval-n,6);
            %Y1NMech is the non-mechanical change in true price (price
            %immediately before execution minus the price immediately after execution of the previous aggressive order
            Y1NMech = AggressiveOrders3(startval+n:endval,6)-AggressiveOrders3(startval:endval-n,13);
            %First with mechanical depletion
            [beta,R,J,~,~,~] = nlinfit(X1,Y1Mech,nlModel,[0;0.15;1]);
            CovMtrxMech = NLLSVarMtrxCalcs(J,R);
            CovMtrxMech = CovMtrxMech/size(R,1);
            betasMech(counter,:) = beta';
            standardErrorsMech(counter,:) = (diag(CovMtrxMech).^0.5)';
            %Now without
            [beta,R,J,~,~,~] = nlinfit(X1,Y1NMech,nlModel,[0;0.15;1]);
            CovMtrxNMech = NLLSVarMtrxCalcs(J,R);
            CovMtrxNMech = CovMtrxNMech/size(R,1);
            betasNMech(counter,:) = beta';
            standardErrorsNMech(counter,:) = (diag(CovMtrxNMech).^0.5)';
            %Now square root estimate
            [betaMech05d,R105d,J05d,~,MSEMech05d,~] = nlinfit(X1,Y1Mech,nlModel2,[0;1]);
            CovMtrxMech05d = NLLSVarMtrxCalcs(J05d,R105d);
            CovMtrxMech05d = CovMtrxMech05d/size(R105d,1);
            betasMech05d(counter,:) = betaMech05d';
            standardErrorsMech05d(counter,:) = (diag(CovMtrxMech05d).^0.5)';
            betaMech05d'
            

            %Tail power law estimation section (optional, not used in paper)
            rets = (AggressiveOrders3(startval+n:endval,6)-AggressiveOrders3(startval:endval-n,6))./AggressiveOrders3(startval:endval-n,6);
            
            

            %Kyle model for robustness (non-mechanical)
            [X,Y,BetaHat,YHat,ehat,Sandwich] = OLSGenericProcessor(X1,Y1NMech);
            z = (BetaHat./diag(sqrt(Sandwich)))*sqrt(size(X1,1));
            r2 = 1- sum(ehat.^2)/sum((Y-mean(Y)).^2);
            lambdasNMech(counter,1) = BetaHat(2,1);
            betasNMechL(counter,:) = BetaHat';
            sandwichesNMech(counter,:) = Sandwich(:)';
            daysizesNMech(counter,1) = size(X1,1);
            r2sNMech(counter,1) = r2;
            
            %Kyle model for robustness (full)
            [X,Y,BetaHat,~,ehat,Sandwich] = OLSGenericProcessor(X1,Y1Mech);
            z = (BetaHat./diag(sqrt(Sandwich)))*sqrt(size(X1,1));
            r2 = 1- sum(ehat.^2)/sum((Y-mean(Y)).^2);
            lambdasMechL(counter,1) = BetaHat(2,1);
            betasMechL(counter,:) = BetaHat';
            sandwichesMechL(counter,:) = Sandwich(:)';
            daysizesMechL(counter,1) = size(X1,1);
            r2sMechL(counter,1) = r2;
            volumesMechL(counter,1) = sum(AggressiveOrders3(startval:endval-n,3));
            additionalValues(counter,:) = [min(AggressiveOrders3(startval:endval-n,6)) max(AggressiveOrders3(startval:endval-n,6)) mean(AggressiveOrders3(startval:endval-n,6)) sum((AggressiveOrders3(startval:endval-n,6)).*(AggressiveOrders3(startval:endval-n,3))) sum(AggressiveOrders3(startval:endval-n,3)) (endval-n+1-startval)];
        
        end
        
        %Store values
        Y1MechAll = [Y1MechAll;AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,6)]; %#ok<AGROW>
        Y1NMechAll = [Y1NMechAll;AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,13)];%#ok<AGROW>
        X1All = [X1All; AggressiveOrders3(1:end-n,3).*AggressiveOrders3(1:end-n,4)];%#ok<AGROW>
        t2All = [t2All; AggressiveOrders3(1:end-n,1)]; %#ok<AGROW>
        retsAll = [retsAll;(AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,6))./AggressiveOrders3(1:end-n,6)]; %#ok<AGROW>
        BetasMech{filecounter,1} = betasMech;
        BetasMech05d{filecounter,1} = betasMech05d;
        StandardErrorsMech{filecounter,1} = standardErrorsMech;
        StandardErrorsMech05d{filecounter,1} = standardErrorsMech05d;
        BetasNMech{filecounter,1} = betasNMech;
        StandardErrorsNMech{filecounter,1} = standardErrorsNMech;


        LambdasNMechL{filecounter,1} = lambdasNMech;
        BetasNMechL{filecounter,1} = betasNMechL;
        SandwichesNMechL{filecounter,1} = sandwichesNMech;
        DaySizesNMechL{filecounter,1} = daysizesNMech;
        R2sNMechL{filecounter,1} = r2sNMech;


        LambdasMechL{filecounter,1} = lambdasMechL;
        BetasMechL{filecounter,1} = betasMechL;
        SandwichesMechL{filecounter,1} = sandwichesMechL;
        DaySizesMechL{filecounter,1} = daysizesMechL;
        R2sMechL{filecounter,1} = r2sMechL;
        VolumesMechL{filecounter,1} = volumesMechL;
        AdditionalValues{filecounter,1} = additionalValues;
        
        DatesForLater{filecounter,1} = datesForLater;
    end

    %Having finished calculating values for all 5 contract months,
    %Cycle through and store values in vectors with daily estimates
    betasExtrasMech = zeros(300,3);
    stErExtrasMech = zeros(300,3);
    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthBetas = size(BetasMech{counter,1},1);
    betasExtrasMech(pointer:pointer+lengthBetas-1,:) = BetasMech{counter,1};
    stErExtrasMech(pointer:pointer+lengthBetas-1,:) = StandardErrorsMech{counter,1};
    pointer = pointer + lengthBetas;
    end
    betasExtrasMech = betasExtrasMech(1:pointer-1,:);
    stErExtrasMech = stErExtrasMech(1:pointer-1,:);
    
    betasExtrasMech05d = zeros(300,2);
    stErExtrasMech05d = zeros(300,2);
    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthBetas = size(BetasMech05d{counter,1},1);
    betasExtrasMech05d(pointer:pointer+lengthBetas-1,:) = BetasMech05d{counter,1};
    stErExtrasMech05d(pointer:pointer+lengthBetas-1,:) = StandardErrorsMech05d{counter,1};
    pointer = pointer + lengthBetas;
    end
    betasExtrasMech05d = betasExtrasMech05d(1:pointer-1,:);
    stErExtrasMech05d = stErExtrasMech05d(1:pointer-1,:);


    betasExtrasNMech = zeros(300,3);
    stErExtrasNMech = zeros(300,3);
    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthBetas = size(BetasNMech{counter,1},1);
    betasExtrasNMech(pointer:pointer+lengthBetas-1,:) = BetasNMech{counter,1};
    stErExtrasNMech(pointer:pointer+lengthBetas-1,:) = StandardErrorsNMech{counter,1};
    pointer = pointer + lengthBetas;
    end
    betasExtrasNMech = betasExtrasNMech(1:pointer-1,:);
    stErExtrasNMech = stErExtrasNMech(1:pointer-1,:);


    lambdaExtrasNoMech = zeros(300,1);
    betaExtrasNoMech = zeros(300,2);
    sandwichesExtrasNoMech = zeros(300,4);
    daySizesNoMech = zeros(300,1);
    r2sNoMech = zeros(300,1);

    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthLambda = size(LambdasNMechL{counter,1},1);
    lambdaExtrasNoMech(pointer:pointer+lengthLambda-1,1) = LambdasNMechL{counter,1};
    betaExtrasNoMech(pointer:pointer+lengthLambda-1,:) = BetasNMechL{counter,1};
    sandwichesExtrasNoMech(pointer:pointer+lengthLambda-1,:) = SandwichesNMechL{counter,1};
    daySizesNoMech(pointer:pointer+lengthLambda-1,:) = DaySizesNMechL{counter,1};
    r2sNoMech(pointer:pointer+lengthLambda-1,:) = R2sNMechL{counter,1};
    pointer = pointer + lengthLambda;
    end
    lambdaExtrasNoMech = lambdaExtrasNoMech(1:pointer-1,:);
    betaExtrasNoMech = betaExtrasNoMech(1:pointer-1,:);
    sandwichesExtrasNoMech = sandwichesExtrasNoMech(1:pointer-1,:);
    daySizesNoMech = daySizesNoMech(1:pointer-1,:);
    r2sNoMech = r2sNoMech(1:pointer-1,:);

    lambdaExtrasMech = zeros(300,1);
    betaExtrasMech = zeros(300,2);
    sandwichesExtrasMech = zeros(300,4);
    daySizesMech = zeros(300,1);
    r2sMech = zeros(300,1);
    volumesMech = zeros(300,1);

    counter = 0;
    pointer = 1;
    while(counter<5)
        counter = counter +1;
        lengthLambda = size(LambdasMechL{counter,1},1);
        lambdaExtrasMech(pointer:pointer+lengthLambda-1,1) = LambdasMechL{counter,1};
        betaExtrasMech(pointer:pointer+lengthLambda-1,:) = BetasMechL{counter,1};
        sandwichesExtrasMech(pointer:pointer+lengthLambda-1,:) = SandwichesMechL{counter,1};
        daySizesMech(pointer:pointer+lengthLambda-1,:) = DaySizesMechL{counter,1};
        r2sMech(pointer:pointer+lengthLambda-1,:) = R2sMechL{counter,1};
        volumesMech(pointer:pointer+lengthLambda-1,:) = VolumesMechL{counter,1};
        pointer = pointer + lengthLambda;
    end

    lambdaExtrasMech = lambdaExtrasMech(1:pointer-1,:);
    betaExtrasMech = betaExtrasMech(1:pointer-1,:);
    sandwichesExtrasMech = sandwichesExtrasMech(1:pointer-1,:);
    daySizesMech = daySizesMech(1:pointer-1,:);
    r2sMech = r2sMech(1:pointer-1,:);
    volumesMech = volumesMech(1:pointer-1,:);

    clearvars a2 AggressiveOrders3 dataIsNaNFormat diffjumps ESU1820180801 RegMtrx2 selector t2H t2M ;
    %plot(betasExtrasMech(:,3))
    

    %Calculate composite non-linear model (full)
    [trimmedData,selector] = VectorTrimmer(hour(t2All),[t2All X1All Y1MechAll Y1NMechAll retsAll],n);
    [betaMechComposite,R1,J,~,MSEMech,~] = nlinfit(trimmedData(:,2),trimmedData(:,3),nlModel,[0;1;1]);
    CovMtrxMech = NLLSVarMtrxCalcs(J,R1);
    CovMtrxMech = CovMtrxMech/size(R1,1);
    betasMech = betaMechComposite';
    standardErrorsMechComposite = (diag(CovMtrxMech).^0.5)';
    betaMechComposite'
    
    %Calculate composite non-linear model (non-mechanical)
    [betaNMechComposite,R2,J,~,MSENMech,~] = nlinfit(trimmedData(:,2),trimmedData(:,4),nlModel,[0;1;1]);
    CovMtrxMech = NLLSVarMtrxCalcs(J,R2);
    CovMtrxMech = CovMtrxMech/size(R2,1);
    (diag(CovMtrxMech).^0.5)';
    betaNMechComposite'
    standardErrorsNMechComposite = (diag(CovMtrxMech).^0.5)';
    (diag(CovMtrxMech).^0.5)'




    
    %Square root estimation section
    [~,Y,BetaHat,YHat,ehat,Sandwich] = OLSGenericProcessor(trimmedData(:,2),trimmedData(:,4));
    zNMechComposite = (BetaHat./diag(sqrt(Sandwich)))*sqrt(size(trimmedData(:,2),1));
    r2 = 1- sum(ehat.^2)/sum((Y-mean(Y)).^2);
    lambdasNMechComposite = BetaHat(2,1);
    betasNMechLComposite = BetaHat';
    sandwichesNMechComposite = Sandwich(:)';
    r2sNMechComposite = r2;
    MSELinearCompositeNMech = mean(ehat.^2);

    [~,Y,BetaHat,~,ehat,Sandwich] = OLSGenericProcessor(trimmedData(:,2),trimmedData(:,3));
    zMechComposite = (BetaHat./diag(sqrt(Sandwich)))*sqrt(size(trimmedData(:,2),1));
    r2 = 1- sum(ehat.^2)/sum((Y-mean(Y)).^2);
    lambdasMechLComposite = BetaHat(2,1);
    betasMechLComposite = BetaHat';
    sandwichesMechLComposite = Sandwich(:)';
    r2sMechLComposite = r2;
    MSELinearCompositeMech = mean(ehat.^2);

    
    syntheticVar = sign(trimmedData(:,2)).*abs(trimmedData(:,2)).^betaNMechComposite(3);
    [~,Y,BetaHat,YHat,ehat,Sandwich] = OLSGenericProcessor(syntheticVar,trimmedData(:,4));
    zNMechComposite = (BetaHat./diag(sqrt(Sandwich)))*sqrt(size(syntheticVar,1));
    r2 = 1- sum(ehat.^2)/sum((Y-mean(Y)).^2);
    lambdasNMechCompositeNL = BetaHat(2,1);
    betasNMechLCompositeNL = BetaHat';
    sandwichesNMechCompositeNL = Sandwich(:)';
    r2sNMechCompositeNL = r2;

    syntheticVar = sign(trimmedData(:,2)).*abs(trimmedData(:,2)).^betaMechComposite(3);
    [~,Y,BetaHat,~,ehat,Sandwich] = OLSGenericProcessor(syntheticVar,trimmedData(:,3));
    zMechComposite = (BetaHat./diag(sqrt(Sandwich)))*sqrt(size(syntheticVar,1));
    r2 = 1- sum(ehat.^2)/sum((Y-mean(Y)).^2);
    lambdasMechLCompositeNL = BetaHat(2,1);
    betasMechLCompositeNL = BetaHat';
    sandwichesMechLCompositeNL = Sandwich(:)';
    r2sMechLCompositeNL = r2;
    

    toc


holidays = [11 28; 1 20; 2 17; 5 31; 7 3; 9 7];
selector = volumesMech>=100000;

%SquareRoot estimation again
nlModel2 = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^0.5).*sign(X))); %Define nonlinear model to estimate
[betaMechComposite05,R105,J05,~,MSEMech05,~] = nlinfit(trimmedData(:,2),trimmedData(:,3),nlModel2,[0;1]);
CovMtrxMech05 = NLLSVarMtrxCalcs(J05,R105);
CovMtrxMech05 = CovMtrxMech05/size(R105,1);
betasMech05 = betaMechComposite05';
standardErrorsMechComposite05 = (diag(CovMtrxMech05).^0.5)';
betaMechComposite05'

%Store values
datesExtras = zeros(300,1);
counter = 0;
pointer = 1;
while(counter<5)
counter = counter +1;
lengthdates = size(DatesForLater{counter,1},1);
datesExtras(pointer:pointer+lengthdates-1,:) = DatesForLater{counter,1};
pointer = pointer + lengthdates;
end
datesExtras = datesExtras(1:pointer-1,:);

%Filter out holidays
selectorNonHolidays = zeros(315,1)+1;
selectorNonHolidays([8 63 127 192 257],:) = 0;
selectorNonHolidays = (selectorNonHolidays>0);
formatOut = 'mm/dd/yy';
dateVec = datestr(datesExtras+1,formatOut);

str = strcat('ComprehensiveStudy_6',Strcatroots{1,outerloop},'All');

%Save to file
save(str);

end

%Redirect back to main directory
cd(parentDirectory);